##Install Packages if Needed
if (!require("ggplot2")) install.packages("ggplot2")
Loading required package: ggplot2
if (!require("effectsize")) install.packages("effectsize")
Loading required package: effectsize
if (!require("emmeans")) install.packages("emmeans")
Loading required package: emmeans
if (!require("dplyr")) install.packages("dplyr")
Loading required package: dplyr
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
if (!require("tidyr")) install.packages("tidyr")
Loading required package: tidyr
if (!require("Rmisc")) install.packages("Rmisc")
Loading required package: Rmisc
Loading required package: lattice
Loading required package: plyr
----------------------------------------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
----------------------------------------------------------------------------------
Attaching package: ‘plyr’
The following objects are masked from ‘package:dplyr’:
arrange, count, desc, failwith, id, mutate, rename, summarise,
summarize
if (!require("ggpubr")) install.packages("ggpubr")
Loading required package: ggpubr
Attaching package: ‘ggpubr’
The following object is masked from ‘package:plyr’:
mutate
if (!require("cowplot")) install.packages("cowplot")
Loading required package: cowplot
Attaching package: ‘cowplot’
The following object is masked from ‘package:ggpubr’:
get_legend
if (!require("gridGraphics")) install.packages("gridGraphics")
Loading required package: gridGraphics
Warning: package ‘gridGraphics’ was built under R version 4.3.2Loading required package: grid
if (!require("tidyr")) install.packages("tidyr")
##Load Packages
library(ggplot2) #Required for plotting
library(effectsize) #Required for eta_squared effect sizes
library(emmeans) #Required for pairwise comparisons
library(dplyr) #Required to seperate columns in dataframe
library(tidyr) #Required for data organization
library(Rmisc) #Required for summarySE for summary statistics
library(ggpubr) #Required for adding pairwise p-values to plots with stat_pvalue_manual
library(cowplot) #Required for arranging ggplots
library(gridGraphics) #Required for adding labels to arranged plots
library(tidyr) #Required for reshaping datafrom from a wide to long format.
Note: Run “Graphing Parameters” section from 01_ExperimentalSetup.R file
Note: Full Data with Bleaching Metrics and Color Scores created in 04_Models.R file
#Load Data
FullData<-read.csv("Outputs/FullData.csv", header=TRUE)
#Set factor variables
FullData$TimeP<-factor(FullData$TimeP, levels=c("W2", "M1", "M4"), ordered=TRUE)
FullData$Site<-factor(FullData$Site, levels=c("KL", "SS"), ordered=TRUE)
FullData$Genotype<-factor(FullData$Genotype, levels=c("AC10", "AC12", "AC8"), ordered=TRUE)
FullData$Treatment<-factor(FullData$Treatment, levels=c("Control", "Heat"), ordered=TRUE)
FullData$Treat<-factor(FullData$Treat, levels=c("C", "H"), ordered=TRUE)
FullData$Seas<-factor(FullData$Seas, levels=c("Summer1", "Summer2", "Winter"), ordered=TRUE)
##Control
FullData_C<-subset(FullData, Treat=="C")
##Heated
FullData_H<-subset(FullData, Treat=="H")
Calculating retention as proportion remaining relative to corresponding control levels (0-1) for each site and genotype at each timepoint.
##Calculate averages for Control Treatment for each Site, Genotype, Timepoint, and Season
names(FullData_C)
[1] "ID" "RandN" "TimeP" "Site" "Genotype"
[6] "Treat" "Treatment" "Seas" "Set" "Score_Full"
[11] "Score_TP" "SA_cm2" "Chl_ug.cm2" "Sym10.6_cm2"
FullData_C.a<-aggregate(FullData_C[,c(10:11, 13:14)], list(FullData_C$Site, FullData_C$Genotype, FullData_C$TimeP, FullData_C$Seas), mean, na.action = na.omit)
names(FullData_C.a)[1:4]<-c("Site", "Genotype", "TimeP", "Seas")
names(FullData_C.a)[5:8]<-paste(names(FullData_C)[c(10:11, 13:14)], "C", sep="_")
##Merge Control Averages with Heated Samples
#Merges by Site, Genotype, Timepoint, and Season
TolData<-merge(FullData_H, FullData_C.a, all.x=TRUE )
##Calculate Proportion Retained for each Bleaching Metric relative to Control
TolData$Score_Full.prop<-round((TolData$Score_Full/TolData$Score_Full_C), 4)
TolData$Score_TP.prop<-round((TolData$Score_TP/TolData$Score_TP_C), 4)
TolData$Chl.prop<-round((TolData$Chl_ug.cm2/TolData$Chl_ug.cm2_C), 4)
TolData$Sym.prop<-round((TolData$Sym10.6_cm2/TolData$Sym10.6_cm2_C), 4)
##Set values >1 to 1
TolData$Score_Full.prop[which(TolData$Score_Full.prop>1)]<-1.0000
TolData$Score_TP.prop[which(TolData$Score_TP.prop>1)]<-1.0000
TolData$Chl.prop[which(TolData$Chl.prop>1)]<-1.0000
TolData$Sym.prop[which(TolData$Sym.prop>1)]<-1.0000
##Create Site and Genotype Variable
TolData$Site.Geno<-paste(TolData$Site, TolData$Genotype, sep="_")
TolData_W2<-subset(TolData, TimeP=="W2")
TolData_M1<-subset(TolData, TimeP=="M1")
TolData_M4<-subset(TolData, TimeP=="M4")
##Check normality
hist(TolData$Chl.prop)
shapiro.test(TolData$Chl.prop)
Shapiro-Wilk normality test
data: TolData$Chl.prop
W = 0.97395, p-value = 0.165
#Normal
##Model as a function of Site and Genotype and Timepoint
Tol_Chl_lm<-lm(Chl.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_lm)); qqline(resid(Tol_Chl_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_lm), resid(Tol_Chl_lm))
#Model Results
summary(Tol_Chl_lm)
Call:
lm(formula = Chl.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.127060 -0.031006 -0.004219 0.036829 0.097690
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.2100753 0.0065828 31.913 < 2e-16 ***
Site.L -0.0093229 0.0092942 -1.003 0.320289
Genotype.L -0.0980285 0.0111500 -8.792 5.31e-12 ***
Genotype.Q -0.0467859 0.0116301 -4.023 0.000180 ***
TimeP.L 0.0638187 0.0112994 5.648 6.22e-07 ***
TimeP.Q 0.0464233 0.0115033 4.036 0.000173 ***
Site.L:Genotype.L 0.0003917 0.0157685 0.025 0.980272
Site.L:Genotype.Q 0.0592529 0.0164183 3.609 0.000673 ***
Site.L:TimeP.L -0.0528050 0.0159798 -3.304 0.001694 **
Site.L:TimeP.Q -0.0236749 0.0161894 -1.462 0.149435
Genotype.L:TimeP.L 0.0406686 0.0194339 2.093 0.041093 *
Genotype.Q:TimeP.L -0.0627603 0.0196595 -3.192 0.002354 **
Genotype.L:TimeP.Q 0.0800244 0.0191901 4.170 0.000111 ***
Genotype.Q:TimeP.Q -0.1328570 0.0206170 -6.444 3.28e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05393 on 54 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.8229, Adjusted R-squared: 0.7803
F-statistic: 19.3 on 13 and 54 DF, p-value: 8.612e-16
anova(Tol_Chl_lm)
Analysis of Variance Table
Response: Chl.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.001563 0.001563 0.5375 0.466633
Genotype 2 0.297307 0.148654 51.1117 3.495e-13 ***
TimeP 2 0.132285 0.066142 22.7418 6.845e-08 ***
Site:Genotype 2 0.041114 0.020557 7.0681 0.001877 **
Site:TimeP 2 0.038643 0.019322 6.6434 0.002634 **
Genotype:TimeP 4 0.218951 0.054738 18.8205 9.691e-10 ***
Residuals 54 0.157054 0.002908
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
----------------------------------------
Site | 1.76e-03 | [0.00, 1.00]
Genotype | 0.34 | [0.16, 1.00]
TimeP | 0.15 | [0.02, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
Site:TimeP | 0.04 | [0.00, 1.00]
Genotype:TimeP | 0.25 | [0.06, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_lm.res<-data.frame(anova(Tol_Chl_lm))
Tol_Chl_lm.res$Predictor<-rownames(Tol_Chl_lm.res)
Tol_Chl_lm.res$EtaSq<-c(eta_squared(Tol_Chl_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_lm.res))
Tol_Chl_lm.res<-Tol_Chl_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with main variables of interest, Site and Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Chl.prop)
shapiro.test(TolData_W2$Chl.prop)
Shapiro-Wilk normality test
data: TolData_W2$Chl.prop
W = 0.98407, p-value = 0.9671
#Normal
##Model as a function of Site and Genotype
Tol_Chl_W2_lm<-lm(Chl.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_W2_lm)); qqline(resid(Tol_Chl_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_W2_lm), resid(Tol_Chl_W2_lm))
#Model Results
summary(Tol_Chl_W2_lm)
Call:
lm(formula = Chl.prop ~ Site + Genotype + Site:Genotype, data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.089900 -0.031587 -0.008183 0.038881 0.103550
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.1842264 0.0119450 15.423 5.03e-11 ***
Site.L 0.0188110 0.0168928 1.114 0.281921
Genotype.L -0.0947759 0.0204291 -4.639 0.000273 ***
Genotype.Q -0.0585870 0.0209464 -2.797 0.012921 *
Site.L:Genotype.L -0.0002833 0.0288911 -0.010 0.992297
Site.L:Genotype.Q 0.0366569 0.0296227 1.237 0.233769
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05552 on 16 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.6758, Adjusted R-squared: 0.5745
F-statistic: 6.671 on 5 and 16 DF, p-value: 0.001551
anova(Tol_Chl_W2_lm)
Analysis of Variance Table
Response: Chl.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.005580 0.005580 1.8106 0.197195
Genotype 2 0.092482 0.046241 15.0038 0.000214 ***
Site:Genotype 2 0.004732 0.002366 0.7677 0.480448
Residuals 16 0.049311 0.003082
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.04 | [0.00, 1.00]
Genotype | 0.61 | [0.29, 1.00]
Site:Genotype | 0.03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_W2_lm.res<-data.frame(anova(Tol_Chl_W2_lm))
Tol_Chl_W2_lm.res$Predictor<-rownames(Tol_Chl_W2_lm.res)
Tol_Chl_W2_lm.res$EtaSq<-c(eta_squared(Tol_Chl_W2_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_W2_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_W2_lm.res))
Tol_Chl_W2_lm.res$TimeP<-rep("W2", nrow(Tol_Chl_W2_lm.res))
Tol_Chl_W2_lm.res<-Tol_Chl_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Chl_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.2033 0.0278 16 0.1445 0.262
AC12 0.2399 0.0278 16 0.1811 0.299
AC8 0.0696 0.0278 16 0.0107 0.128
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.2514 0.0278 16 0.1925 0.310
AC12 0.2242 0.0321 16 0.1563 0.292
AC8 0.1170 0.0321 16 0.0491 0.185
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.0366 0.0393 16 -0.933 0.6280
AC10 - AC8 0.1338 0.0393 16 3.407 0.0095
AC12 - AC8 0.1704 0.0393 16 4.340 0.0014
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0272 0.0424 16 0.640 0.8003
AC10 - AC8 0.1343 0.0424 16 3.168 0.0156
AC12 - AC8 0.1072 0.0453 16 2.364 0.0753
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Chl_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.2033 0.0278 16 0.1445 0.262
SS 0.2514 0.0278 16 0.1925 0.310
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.2399 0.0278 16 0.1811 0.299
SS 0.2242 0.0321 16 0.1563 0.292
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.0696 0.0278 16 0.0107 0.128
SS 0.1170 0.0321 16 0.0491 0.185
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0481 0.0393 16 -1.224 0.2387
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.0157 0.0424 16 0.371 0.7156
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.0475 0.0424 16 -1.120 0.2793
##Save p-values
#Genotypes within Sites
Tol_Chl_W2_lm.geno<-data.frame(emmeans(Tol_Chl_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_W2_lm.geno<-Tol_Chl_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_W2_lm.geno$group1<-paste(Tol_Chl_W2_lm.geno$Site, Tol_Chl_W2_lm.geno$group1, sep="_")
Tol_Chl_W2_lm.geno$group2<-paste(Tol_Chl_W2_lm.geno$Site, Tol_Chl_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Chl_W2_lm.site<-data.frame(emmeans(Tol_Chl_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_W2_lm.site<-Tol_Chl_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_W2_lm.site$group1<-paste(Tol_Chl_W2_lm.site$group1, Tol_Chl_W2_lm.site$Genotype, sep="_")
Tol_Chl_W2_lm.site$group2<-paste(Tol_Chl_W2_lm.site$group2, Tol_Chl_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Chl_W2_lm.p<-rbind(Tol_Chl_W2_lm.geno[,c(1:2,4:8)], Tol_Chl_W2_lm.site[,c(1:2,4:8)])
Tol_Chl_W2_lm.p<-Tol_Chl_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Chl_W2_lm.p$Sig<-ifelse(Tol_Chl_W2_lm.p$p<0.001, "***", ifelse(Tol_Chl_W2_lm.p$p<0.01, "**", ifelse(Tol_Chl_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Chl_W2_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_W2_lm.p))
Tol_Chl_W2_lm.p$TimeP<-rep("W2", nrow(Tol_Chl_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_Chl_W2_SG<-summarySE(TolData_W2, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Chl_W2_SG.plot<-ggplot(Tol_Chl_W2_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Chlorophyll Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_W2_lm.p, y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE); Tol_Chl_W2_SG.plot
##Check normality
hist(TolData_M1$Chl.prop)
shapiro.test(TolData_M1$Chl.prop)
Shapiro-Wilk normality test
data: TolData_M1$Chl.prop
W = 0.9306, p-value = 0.1263
#Normal
##Model as a function of Site and Genotype
Tol_Chl_M1_lm<-lm(Chl.prop~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_M1_lm)); qqline(resid(Tol_Chl_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_M1_lm), resid(Tol_Chl_M1_lm))
#Model Results
summary(Tol_Chl_M1_lm)
Call:
lm(formula = Chl.prop ~ Site + Genotype + Site:Genotype, data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.107050 -0.034575 0.003479 0.025531 0.088050
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.17217 0.01178 14.619 1.12e-10 ***
Site.L 0.01036 0.01666 0.622 0.5428
Genotype.L -0.16337 0.01935 -8.442 2.74e-07 ***
Genotype.Q 0.06169 0.02139 2.883 0.0108 *
Site.L:Genotype.L 0.03819 0.02737 1.395 0.1820
Site.L:Genotype.Q 0.05454 0.03026 1.803 0.0903 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05474 on 16 degrees of freedom
Multiple R-squared: 0.8424, Adjusted R-squared: 0.7932
F-statistic: 17.11 on 5 and 16 DF, p-value: 6.411e-06
anova(Tol_Chl_M1_lm)
Analysis of Variance Table
Response: Chl.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.002283 0.002283 0.7619 0.3956
Genotype 2 0.238424 0.119212 39.7888 6.168e-07 ***
Site:Genotype 2 0.015569 0.007785 2.5982 0.1054
Residuals 16 0.047938 0.002996
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 7.50e-03 | [0.00, 1.00]
Genotype | 0.78 | [0.58, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_M1_lm.res<-data.frame(anova(Tol_Chl_M1_lm))
Tol_Chl_M1_lm.res$Predictor<-rownames(Tol_Chl_M1_lm.res)
Tol_Chl_M1_lm.res$EtaSq<-c(eta_squared(Tol_Chl_M1_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_M1_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_M1_lm.res))
Tol_Chl_M1_lm.res$TimeP<-rep("M1", nrow(Tol_Chl_M1_lm.res))
Tol_Chl_M1_lm.res<-Tol_Chl_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Chl_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.3089 0.0274 16 0.2509 0.3669
AC12 0.1460 0.0316 16 0.0790 0.2130
AC8 0.0397 0.0274 16 -0.0183 0.0977
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.3169 0.0274 16 0.2588 0.3749
AC12 0.0976 0.0316 16 0.0306 0.1646
AC8 0.1240 0.0274 16 0.0660 0.1820
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1629 0.0418 16 3.897 0.0035
AC10 - AC8 0.2692 0.0387 16 6.956 <.0001
AC12 - AC8 0.1063 0.0418 16 2.543 0.0538
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.2192 0.0418 16 5.244 0.0002
AC10 - AC8 0.1928 0.0387 16 4.983 0.0004
AC12 - AC8 -0.0264 0.0418 16 -0.631 0.8056
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Chl_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.3089 0.0274 16 0.2509 0.3669
SS 0.3169 0.0274 16 0.2588 0.3749
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.1460 0.0316 16 0.0790 0.2130
SS 0.0976 0.0316 16 0.0306 0.1646
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.0397 0.0274 16 -0.0183 0.0977
SS 0.1240 0.0274 16 0.0660 0.1820
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.00795 0.0387 16 -0.205 0.8398
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.04833 0.0447 16 1.081 0.2955
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.08432 0.0387 16 -2.179 0.0447
##Save p-values
#Genotypes within Sites
Tol_Chl_M1_lm.geno<-data.frame(emmeans(Tol_Chl_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_M1_lm.geno<-Tol_Chl_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M1_lm.geno$group1<-paste(Tol_Chl_M1_lm.geno$Site, Tol_Chl_M1_lm.geno$group1, sep="_")
Tol_Chl_M1_lm.geno$group2<-paste(Tol_Chl_M1_lm.geno$Site, Tol_Chl_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Chl_M1_lm.site<-data.frame(emmeans(Tol_Chl_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_M1_lm.site<-Tol_Chl_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M1_lm.site$group1<-paste(Tol_Chl_M1_lm.site$group1, Tol_Chl_M1_lm.site$Genotype, sep="_")
Tol_Chl_M1_lm.site$group2<-paste(Tol_Chl_M1_lm.site$group2, Tol_Chl_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Chl_M1_lm.p<-rbind(Tol_Chl_M1_lm.geno[,c(1:2,4:8)], Tol_Chl_M1_lm.site[,c(1:2,4:8)])
Tol_Chl_M1_lm.p<-Tol_Chl_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Chl_M1_lm.p$Sig<-ifelse(Tol_Chl_M1_lm.p$p<0.001, "***", ifelse(Tol_Chl_M1_lm.p$p<0.01, "**", ifelse(Tol_Chl_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Chl_M1_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_M1_lm.p))
Tol_Chl_M1_lm.p$TimeP<-rep("M1", nrow(Tol_Chl_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_Chl_M1_SG<-summarySE(TolData_M1, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Chl_M1_SG.plot<-ggplot(Tol_Chl_M1_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Chlorophyll Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_M1_lm.p, y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE); Tol_Chl_M1_SG.plot
##Check normality
hist(TolData_M4$Chl.prop)
shapiro.test(TolData_M4$Chl.prop)
Shapiro-Wilk normality test
data: TolData_M4$Chl.prop
W = 0.91299, p-value = 0.04096
#Not Normal
##Try log+1 transformation
hist(log(TolData_M4$Chl.prop+1))
shapiro.test(log(TolData_M4$Chl.prop+1))
Shapiro-Wilk normality test
data: log(TolData_M4$Chl.prop + 1)
W = 0.93102, p-value = 0.1028
#Normal
##Model as a function of Site and Genotype
##Model with log transformation
Tol_Chl_M4_lm<-lm(log(Chl.prop+1)~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_M4_lm)); qqline(resid(Tol_Chl_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_M4_lm), resid(Tol_Chl_M4_lm))
#Model Results
summary(Tol_Chl_M4_lm)
Call:
lm(formula = log(Chl.prop + 1) ~ Site + Genotype + Site:Genotype,
data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.062753 -0.025900 -0.004896 0.017009 0.065759
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.23863 0.00807 29.568 < 2e-16 ***
Site.L -0.04272 0.01141 -3.743 0.00149 **
Genotype.L -0.03027 0.01398 -2.166 0.04400 *
Genotype.Q -0.11052 0.01398 -7.907 2.9e-07 ***
Site.L:Genotype.L -0.03189 0.01977 -1.613 0.12415
Site.L:Genotype.Q 0.05753 0.01977 2.910 0.00934 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03954 on 18 degrees of freedom
Multiple R-squared: 0.8368, Adjusted R-squared: 0.7915
F-statistic: 18.46 on 5 and 18 DF, p-value: 1.598e-06
anova(Tol_Chl_M4_lm)
Analysis of Variance Table
Response: log(Chl.prop + 1)
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.021901 0.021901 14.0102 0.001489 **
Genotype 2 0.105056 0.052528 33.6031 8.379e-07 ***
Site:Genotype 2 0.017304 0.008652 5.5349 0.013381 *
Residuals 18 0.028137 0.001563
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.13 | [0.00, 1.00]
Genotype | 0.61 | [0.32, 1.00]
Site:Genotype | 0.10 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_M4_lm.res<-data.frame(anova(Tol_Chl_M4_lm))
Tol_Chl_M4_lm.res$Predictor<-rownames(Tol_Chl_M4_lm.res)
Tol_Chl_M4_lm.res$EtaSq<-c(eta_squared(Tol_Chl_M4_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_M4_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_M4_lm.res))
Tol_Chl_M4_lm.res$TimeP<-rep("M4", nrow(Tol_Chl_M4_lm.res))
Tol_Chl_M4_lm.res<-Tol_Chl_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype and Site*Genotype.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Chl_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.213 0.0198 18 0.171 0.254
AC12 0.392 0.0198 18 0.351 0.434
AC8 0.202 0.0198 18 0.160 0.243
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.217 0.0198 18 0.176 0.259
AC12 0.265 0.0198 18 0.224 0.307
AC8 0.143 0.0198 18 0.101 0.184
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.1797 0.028 18 -6.428 <.0001
AC10 - AC8 0.0109 0.028 18 0.391 0.9196
AC12 - AC8 0.1906 0.028 18 6.819 <.0001
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.0482 0.028 18 -1.724 0.2237
AC10 - AC8 0.0747 0.028 18 2.672 0.0393
AC12 - AC8 0.1229 0.028 18 4.396 0.0010
Note: contrasts are still on the log(mu + 1) scale
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Chl_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.213 0.0198 18 0.171 0.254
SS 0.217 0.0198 18 0.176 0.259
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.392 0.0198 18 0.351 0.434
SS 0.265 0.0198 18 0.224 0.307
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.202 0.0198 18 0.160 0.243
SS 0.143 0.0198 18 0.101 0.184
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.00468 0.028 18 -0.167 0.8688
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.12684 0.028 18 4.537 0.0003
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.05909 0.028 18 2.114 0.0488
Note: contrasts are still on the log(mu + 1) scale
##Save p-values
#Genotypes within Sites
Tol_Chl_M4_lm.geno<-data.frame(emmeans(Tol_Chl_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_M4_lm.geno<-Tol_Chl_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M4_lm.geno$group1<-paste(Tol_Chl_M4_lm.geno$Site, Tol_Chl_M4_lm.geno$group1, sep="_")
Tol_Chl_M4_lm.geno$group2<-paste(Tol_Chl_M4_lm.geno$Site, Tol_Chl_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Chl_M4_lm.site<-data.frame(emmeans(Tol_Chl_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_M4_lm.site<-Tol_Chl_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M4_lm.site$group1<-paste(Tol_Chl_M4_lm.site$group1, Tol_Chl_M4_lm.site$Genotype, sep="_")
Tol_Chl_M4_lm.site$group2<-paste(Tol_Chl_M4_lm.site$group2, Tol_Chl_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Chl_M4_lm.p<-rbind(Tol_Chl_M4_lm.geno[,c(1:2,4:8)], Tol_Chl_M4_lm.site[,c(1:2,4:8)])
Tol_Chl_M4_lm.p<-Tol_Chl_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Chl_M4_lm.p$Sig<-ifelse(Tol_Chl_M4_lm.p$p<0.001, "***", ifelse(Tol_Chl_M4_lm.p$p<0.01, "**", ifelse(Tol_Chl_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Chl_M4_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_M4_lm.p))
Tol_Chl_M4_lm.p$TimeP<-rep("M4", nrow(Tol_Chl_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_Chl_M4_SG<-summarySE(TolData_M4, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Chl_M4_SG.plot<-ggplot(Tol_Chl_M4_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Chlorophyll Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_M4_lm.p, y.position=0.55, step.increase=0.20, label="Sig", hide.ns=TRUE); Tol_Chl_M4_SG.plot
##Check normality
hist(TolData$Sym.prop)
shapiro.test(TolData$Sym.prop)
Shapiro-Wilk normality test
data: TolData$Sym.prop
W = 0.92235, p-value = 0.0003707
#Not Normal
##Try log+1 transformation
hist(log(TolData$Sym.prop+1))
shapiro.test(log(TolData$Sym.prop+1))
Shapiro-Wilk normality test
data: log(TolData$Sym.prop + 1)
W = 0.92164, p-value = 0.000345
#Not normal
##Try square transformation
hist((TolData$Sym.prop)^2)
shapiro.test((TolData$Sym.prop)^2)
Shapiro-Wilk normality test
data: (TolData$Sym.prop)^2
W = 0.87147, p-value = 3.933e-06
#Not normal
##Try cubed transformation
hist((TolData$Sym.prop)^3)
shapiro.test((TolData$Sym.prop)^3)
Shapiro-Wilk normality test
data: (TolData$Sym.prop)^3
W = 0.80697, p-value = 4.342e-08
#Not normal
##Model as a function of Site and Genotype and Timepoint
##Model with no transformation and check residuals
Tol_Sym_lm<-lm(Sym.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_lm)); qqline(resid(Tol_Sym_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_lm), resid(Tol_Sym_lm))
#Model Results
summary(Tol_Sym_lm)
Call:
lm(formula = Sym.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.37464 -0.08630 -0.00192 0.08972 0.35883
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.544513 0.018576 29.313 < 2e-16 ***
Site.L -0.001834 0.026197 -0.070 0.944443
Genotype.L -0.373179 0.031737 -11.758 < 2e-16 ***
Genotype.Q 0.028101 0.032605 0.862 0.392497
TimeP.L 0.133503 0.031737 4.206 9.64e-05 ***
TimeP.Q 0.190344 0.032605 5.838 2.94e-07 ***
Site.L:Genotype.L -0.060389 0.044883 -1.345 0.183999
Site.L:Genotype.Q 0.074992 0.045930 1.633 0.108240
Site.L:TimeP.L -0.159707 0.044883 -3.558 0.000778 ***
Site.L:TimeP.Q -0.027928 0.045930 -0.608 0.545653
Genotype.L:TimeP.L 0.085371 0.055316 1.543 0.128489
Genotype.Q:TimeP.L -0.084924 0.054623 -1.555 0.125749
Genotype.L:TimeP.Q 0.065986 0.054623 1.208 0.232204
Genotype.Q:TimeP.Q -0.205029 0.058264 -3.519 0.000878 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1535 on 55 degrees of freedom
Multiple R-squared: 0.8036, Adjusted R-squared: 0.7572
F-statistic: 17.31 on 13 and 55 DF, p-value: 6.022e-15
anova(Tol_Sym_lm)
Analysis of Variance Table
Response: Sym.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.0007 0.00068 0.0289 0.865708
Genotype 2 3.2796 1.63981 69.5880 8.602e-16 ***
TimeP 2 1.1556 0.57780 24.5197 2.438e-08 ***
Site:Genotype 2 0.1111 0.05557 2.3581 0.104096
Site:TimeP 2 0.3138 0.15688 6.6576 0.002575 **
Genotype:TimeP 4 0.4428 0.11069 4.6974 0.002474 **
Residuals 55 1.2961 0.02356
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
----------------------------------------
Site | 1.03e-04 | [0.00, 1.00]
Genotype | 0.50 | [0.33, 1.00]
TimeP | 0.18 | [0.04, 1.00]
Site:Genotype | 0.02 | [0.00, 1.00]
Site:TimeP | 0.05 | [0.00, 1.00]
Genotype:TimeP | 0.07 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_lm.res<-data.frame(anova(Tol_Sym_lm))
Tol_Sym_lm.res$Predictor<-rownames(Tol_Sym_lm.res)
Tol_Sym_lm.res$EtaSq<-c(eta_squared(Tol_Sym_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_lm.res))
Tol_Sym_lm.res<-Tol_Sym_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with main variables of interest, Site and Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Sym.prop)
shapiro.test(TolData_W2$Sym.prop)
Shapiro-Wilk normality test
data: TolData_W2$Sym.prop
W = 0.93605, p-value = 0.1477
#Normal
##Model as a function of Site and Genotype
Tol_Sym_W2_lm<-lm(Sym.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_W2_lm)); qqline(resid(Tol_Sym_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_W2_lm), resid(Tol_Sym_W2_lm))
#Model Results
summary(Tol_Sym_W2_lm)
Call:
lm(formula = Sym.prop ~ Site + Genotype + Site:Genotype, data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.26060 -0.11145 0.00000 0.08525 0.28580
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5237625 0.0325072 16.112 9.90e-12 ***
Site.L 0.0939568 0.0459722 2.044 0.0568 .
Genotype.L -0.4152131 0.0570402 -7.279 1.29e-06 ***
Genotype.Q -0.0005205 0.0555584 -0.009 0.9926
Site.L:Genotype.L -0.2061000 0.0806671 -2.555 0.0205 *
Site.L:Genotype.Q 0.0462891 0.0785715 0.589 0.5635
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.155 on 17 degrees of freedom
Multiple R-squared: 0.7924, Adjusted R-squared: 0.7313
F-statistic: 12.98 on 5 and 17 DF, p-value: 2.647e-05
anova(Tol_Sym_W2_lm)
Analysis of Variance Table
Response: Sym.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.16148 0.16148 6.7208 0.01897 *
Genotype 2 1.22843 0.61421 25.5640 7.508e-06 ***
Site:Genotype 2 0.16883 0.08441 3.5133 0.05283 .
Residuals 17 0.40845 0.02403
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.08 | [0.00, 1.00]
Genotype | 0.62 | [0.33, 1.00]
Site:Genotype | 0.09 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_W2_lm.res<-data.frame(anova(Tol_Sym_W2_lm))
Tol_Sym_W2_lm.res$Predictor<-rownames(Tol_Sym_W2_lm.res)
Tol_Sym_W2_lm.res$EtaSq<-c(eta_squared(Tol_Sym_W2_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_W2_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_W2_lm.res))
Tol_Sym_W2_lm.res$TimeP<-rep("W2", nrow(Tol_Sym_W2_lm.res))
Tol_Sym_W2_lm.res<-Tol_Sym_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Marginal effect (p<0.1) of Site * Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Sym_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.634 0.0775 17 0.4708 0.798
AC12 0.484 0.0775 17 0.3210 0.648
AC8 0.253 0.0775 17 0.0897 0.417
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 1.000 0.0775 17 0.8365 1.164
AC12 0.564 0.0775 17 0.4004 0.727
AC8 0.207 0.0895 17 0.0179 0.396
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.150 0.110 17 1.367 0.3795
AC10 - AC8 0.381 0.110 17 3.477 0.0077
AC12 - AC8 0.231 0.110 17 2.110 0.1174
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.436 0.110 17 3.979 0.0026
AC10 - AC8 0.793 0.118 17 6.701 <.0001
AC12 - AC8 0.357 0.118 17 3.017 0.0201
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Sym_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.634 0.0775 17 0.4708 0.798
SS 1.000 0.0775 17 0.8365 1.164
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.484 0.0775 17 0.3210 0.648
SS 0.564 0.0775 17 0.4004 0.727
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.253 0.0775 17 0.0897 0.417
SS 0.207 0.0895 17 0.0179 0.396
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.3657 0.110 17 -3.337 0.0039
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0794 0.110 17 -0.725 0.4785
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.0465 0.118 17 0.393 0.6994
##Save p-values
#Genotypes within Sites
Tol_Sym_W2_lm.geno<-data.frame(emmeans(Tol_Sym_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_W2_lm.geno<-Tol_Sym_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_W2_lm.geno$group1<-paste(Tol_Sym_W2_lm.geno$Site, Tol_Sym_W2_lm.geno$group1, sep="_")
Tol_Sym_W2_lm.geno$group2<-paste(Tol_Sym_W2_lm.geno$Site, Tol_Sym_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Sym_W2_lm.site<-data.frame(emmeans(Tol_Sym_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_W2_lm.site<-Tol_Sym_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_W2_lm.site$group1<-paste(Tol_Sym_W2_lm.site$group1, Tol_Sym_W2_lm.site$Genotype, sep="_")
Tol_Sym_W2_lm.site$group2<-paste(Tol_Sym_W2_lm.site$group2, Tol_Sym_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Sym_W2_lm.p<-rbind(Tol_Sym_W2_lm.geno[,c(1:2,4:8)], Tol_Sym_W2_lm.site[,c(1:2,4:8)])
Tol_Sym_W2_lm.p<-Tol_Sym_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Sym_W2_lm.p$Sig<-ifelse(Tol_Sym_W2_lm.p$p<0.001, "***", ifelse(Tol_Sym_W2_lm.p$p<0.01, "**", ifelse(Tol_Sym_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Sym_W2_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_W2_lm.p))
Tol_Sym_W2_lm.p$TimeP<-rep("W2", nrow(Tol_Sym_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_Sym_W2_SG<-summarySE(TolData_W2, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Sym_W2_SG.plot<-ggplot(Tol_Sym_W2_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Symbiont Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_W2_lm.p, y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_W2_SG.plot
##Check normality
hist(TolData_M1$Sym.prop)
shapiro.test(TolData_M1$Sym.prop)
Shapiro-Wilk normality test
data: TolData_M1$Sym.prop
W = 0.85748, p-value = 0.004618
#Not normal
##Try log+1 transformation
hist(log(TolData_M1$Sym.prop+1))
shapiro.test(log(TolData_M1$Sym.prop+1))
Shapiro-Wilk normality test
data: log(TolData_M1$Sym.prop + 1)
W = 0.87335, p-value = 0.009041
#Not normal but improved
##Model as a function of Site and Genotype
##Model with log transformation
Tol_Sym_M1_lm<-lm(log(Sym.prop+1)~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_M1_lm)); qqline(resid(Tol_Sym_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_M1_lm), resid(Tol_Sym_M1_lm))
#Model Results
summary(Tol_Sym_M1_lm)
Call:
lm(formula = log(Sym.prop + 1) ~ Site + Genotype + Site:Genotype,
data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.098079 -0.055799 -0.004888 0.043331 0.151113
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.307716 0.016556 18.587 2.95e-12 ***
Site.L 0.027767 0.023413 1.186 0.252954
Genotype.L -0.295494 0.027204 -10.862 8.59e-09 ***
Genotype.Q 0.129387 0.030075 4.302 0.000548 ***
Site.L:Genotype.L 0.080682 0.038472 2.097 0.052224 .
Site.L:Genotype.Q -0.005572 0.042532 -0.131 0.897412
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.07694 on 16 degrees of freedom
Multiple R-squared: 0.8989, Adjusted R-squared: 0.8673
F-statistic: 28.46 on 5 and 16 DF, p-value: 2e-07
anova(Tol_Sym_M1_lm)
Analysis of Variance Table
Response: log(Sym.prop + 1)
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.00823 0.00823 1.3902 0.2556
Genotype 2 0.80811 0.40406 68.2486 1.468e-08 ***
Site:Genotype 2 0.02614 0.01307 2.2076 0.1423
Residuals 16 0.09473 0.00592
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 8.78e-03 | [0.00, 1.00]
Genotype | 0.86 | [0.73, 1.00]
Site:Genotype | 0.03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_M1_lm.res<-data.frame(anova(Tol_Sym_M1_lm))
Tol_Sym_M1_lm.res$Predictor<-rownames(Tol_Sym_M1_lm.res)
Tol_Sym_M1_lm.res$EtaSq<-c(eta_squared(Tol_Sym_M1_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_M1_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_M1_lm.res))
Tol_Sym_M1_lm.res$TimeP<-rep("M1", nrow(Tol_Sym_M1_lm.res))
Tol_Sym_M1_lm.res<-Tol_Sym_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Sym_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.5918 0.0385 16 0.5102 0.673
AC12 0.1792 0.0444 16 0.0850 0.273
AC8 0.0932 0.0385 16 0.0117 0.175
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.5472 0.0385 16 0.4656 0.629
AC12 0.2249 0.0444 16 0.1307 0.319
AC8 0.2100 0.0385 16 0.1284 0.292
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.413 0.0588 16 7.021 <.0001
AC10 - AC8 0.499 0.0544 16 9.164 <.0001
AC12 - AC8 0.086 0.0588 16 1.463 0.3339
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.322 0.0588 16 5.483 0.0001
AC10 - AC8 0.337 0.0544 16 6.198 <.0001
AC12 - AC8 0.015 0.0588 16 0.255 0.9650
Note: contrasts are still on the log(mu + 1) scale
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Sym_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.5918 0.0385 16 0.5102 0.673
SS 0.5472 0.0385 16 0.4656 0.629
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.1792 0.0444 16 0.0850 0.273
SS 0.2249 0.0444 16 0.1307 0.319
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.0932 0.0385 16 0.0117 0.175
SS 0.2100 0.0385 16 0.1284 0.292
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS 0.0446 0.0544 16 0.820 0.4241
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0457 0.0628 16 -0.727 0.4775
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.1167 0.0544 16 -2.146 0.0476
Note: contrasts are still on the log(mu + 1) scale
##Save p-values
#Genotypes within Sites
Tol_Sym_M1_lm.geno<-data.frame(emmeans(Tol_Sym_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_M1_lm.geno<-Tol_Sym_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M1_lm.geno$group1<-paste(Tol_Sym_M1_lm.geno$Site, Tol_Sym_M1_lm.geno$group1, sep="_")
Tol_Sym_M1_lm.geno$group2<-paste(Tol_Sym_M1_lm.geno$Site, Tol_Sym_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Sym_M1_lm.site<-data.frame(emmeans(Tol_Sym_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_M1_lm.site<-Tol_Sym_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M1_lm.site$group1<-paste(Tol_Sym_M1_lm.site$group1, Tol_Sym_M1_lm.site$Genotype, sep="_")
Tol_Sym_M1_lm.site$group2<-paste(Tol_Sym_M1_lm.site$group2, Tol_Sym_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Sym_M1_lm.p<-rbind(Tol_Sym_M1_lm.geno[,c(1:2,4:8)], Tol_Sym_M1_lm.site[,c(1:2,4:8)])
Tol_Sym_M1_lm.p<-Tol_Sym_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Sym_M1_lm.p$Sig<-ifelse(Tol_Sym_M1_lm.p$p<0.001, "***", ifelse(Tol_Sym_M1_lm.p$p<0.01, "**", ifelse(Tol_Sym_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Sym_M1_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_M1_lm.p))
Tol_Sym_M1_lm.p$TimeP<-rep("M1", nrow(Tol_Sym_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_Sym_M1_SG<-summarySE(TolData_M1, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Sym_M1_SG.plot<-ggplot(Tol_Sym_M1_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Symbiont Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_M1_lm.p, y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_M1_SG.plot
##Check normality
hist(TolData_M4$Sym.prop)
shapiro.test(TolData_M4$Sym.prop)
Shapiro-Wilk normality test
data: TolData_M4$Sym.prop
W = 0.88557, p-value = 0.01078
#Not Normal
##Try square transformation
hist((TolData_M4$Sym.prop)^2)
shapiro.test((TolData_M4$Sym.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M4$Sym.prop)^2
W = 0.87973, p-value = 0.008208
#Not Normal
##Try cubed transformation
hist((TolData_M4$Sym.prop)^3)
shapiro.test((TolData_M4$Sym.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M4$Sym.prop)^3
W = 0.8488, p-value = 0.002077
#Not Normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_Sym_M4_lm<-lm(Sym.prop~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_M4_lm)); qqline(resid(Tol_Sym_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_M4_lm), resid(Tol_Sym_M4_lm))
#Model Results
summary(Tol_Sym_M4_lm)
Call:
lm(formula = Sym.prop ~ Site + Genotype + Site:Genotype, data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.31820 -0.08804 0.00000 0.08618 0.29740
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.71662 0.03157 22.702 1.07e-14 ***
Site.L -0.12617 0.04464 -2.826 0.0112 *
Genotype.L -0.28587 0.05467 -5.229 5.67e-05 ***
Genotype.Q -0.11565 0.05467 -2.115 0.0486 *
Site.L:Genotype.L -0.09494 0.07732 -1.228 0.2353
Site.L:Genotype.Q 0.17275 0.07732 2.234 0.0384 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1546 on 18 degrees of freedom
Multiple R-squared: 0.7201, Adjusted R-squared: 0.6423
F-statistic: 9.26 on 5 and 18 DF, p-value: 0.0001686
anova(Tol_Sym_M4_lm)
Analysis of Variance Table
Response: Sym.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.19101 0.19101 7.9874 0.011190 *
Genotype 2 0.76080 0.38040 15.9068 0.000105 ***
Site:Genotype 2 0.15542 0.07771 3.2496 0.062386 .
Residuals 18 0.43046 0.02391
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.12 | [0.00, 1.00]
Genotype | 0.49 | [0.18, 1.00]
Site:Genotype | 0.10 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_M4_lm.res<-data.frame(anova(Tol_Sym_M4_lm))
Tol_Sym_M4_lm.res$Predictor<-rownames(Tol_Sym_M4_lm.res)
Tol_Sym_M4_lm.res$EtaSq<-c(eta_squared(Tol_Sym_M4_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_M4_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_M4_lm.res))
Tol_Sym_M4_lm.res$TimeP<-rep("M4", nrow(Tol_Sym_M4_lm.res))
Tol_Sym_M4_lm.res<-Tol_Sym_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Marginal effect (p<0.1) of Site * Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Sym_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.863 0.0773 18 0.701 1.026
AC12 1.000 0.0773 18 0.838 1.162
AC8 0.554 0.0773 18 0.392 0.717
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.880 0.0773 18 0.717 1.042
AC12 0.622 0.0773 18 0.460 0.785
AC8 0.380 0.0773 18 0.218 0.543
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.137 0.109 18 -1.249 0.4409
AC10 - AC8 0.309 0.109 18 2.829 0.0285
AC12 - AC8 0.446 0.109 18 4.078 0.0019
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.258 0.109 18 2.356 0.0734
AC10 - AC8 0.499 0.109 18 4.565 0.0007
AC12 - AC8 0.242 0.109 18 2.210 0.0965
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Sym_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.863 0.0773 18 0.701 1.026
SS 0.880 0.0773 18 0.717 1.042
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 1.000 0.0773 18 0.838 1.162
SS 0.622 0.0773 18 0.460 0.785
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.554 0.0773 18 0.392 0.717
SS 0.380 0.0773 18 0.218 0.543
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0163 0.109 18 -0.149 0.8835
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.3779 0.109 18 3.456 0.0028
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.1736 0.109 18 1.588 0.1297
##Save p-values
#Genotypes within Sites
Tol_Sym_M4_lm.geno<-data.frame(emmeans(Tol_Sym_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_M4_lm.geno<-Tol_Sym_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M4_lm.geno$group1<-paste(Tol_Sym_M4_lm.geno$Site, Tol_Sym_M4_lm.geno$group1, sep="_")
Tol_Sym_M4_lm.geno$group2<-paste(Tol_Sym_M4_lm.geno$Site, Tol_Sym_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Sym_M4_lm.site<-data.frame(emmeans(Tol_Sym_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_M4_lm.site<-Tol_Sym_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M4_lm.site$group1<-paste(Tol_Sym_M4_lm.site$group1, Tol_Sym_M4_lm.site$Genotype, sep="_")
Tol_Sym_M4_lm.site$group2<-paste(Tol_Sym_M4_lm.site$group2, Tol_Sym_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Sym_M4_lm.p<-rbind(Tol_Sym_M4_lm.geno[,c(1:2,4:8)], Tol_Sym_M4_lm.site[,c(1:2,4:8)])
Tol_Sym_M4_lm.p<-Tol_Sym_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Sym_M4_lm.p$Sig<-ifelse(Tol_Sym_M4_lm.p$p<0.001, "***", ifelse(Tol_Sym_M4_lm.p$p<0.01, "**", ifelse(Tol_Sym_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Sym_M4_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_M4_lm.p))
Tol_Sym_M4_lm.p$TimeP<-rep("M4", nrow(Tol_Sym_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_Sym_M4_SG<-summarySE(TolData_M4, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Sym_M4_SG.plot<-ggplot(Tol_Sym_M4_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Symbiont Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_M4_lm.p, y.position=1.1, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_M4_SG.plot
##Check normality
hist(TolData$Score_Full.prop)
shapiro.test(TolData$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData$Score_Full.prop
W = 0.90488, p-value = 6.821e-05
#Not Normal
##Try square transformation
hist((TolData$Score_Full.prop)^2)
shapiro.test((TolData$Score_Full.prop)^2)
Shapiro-Wilk normality test
data: (TolData$Score_Full.prop)^2
W = 0.9194, p-value = 0.0002754
#Not normal
##Try cubed transformation
hist((TolData$Score_Full.prop)^3)
shapiro.test((TolData$Score_Full.prop)^3)
Shapiro-Wilk normality test
data: (TolData$Score_Full.prop)^3
W = 0.91645, p-value = 0.0002055
#Not normal
##Model as a function of Site and Genotype and Timepoint
##Model with no transformation and check residuals
Tol_ScoreF_lm<-lm(Score_Full.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_lm)); qqline(resid(Tol_ScoreF_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_lm), resid(Tol_ScoreF_lm))
Residuals are not great, need to check for other modeling options for Color Full Set.
#Model Results
summary(Tol_ScoreF_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.217253 -0.020809 0.001419 0.037990 0.181147
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.816721 0.010646 76.716 < 2e-16 ***
Site.L 0.068538 0.015014 4.565 2.85e-05 ***
Genotype.L -0.162107 0.018189 -8.912 2.92e-12 ***
Genotype.Q -0.000730 0.018686 -0.039 0.96898
TimeP.L 0.143092 0.018189 7.867 1.44e-10 ***
TimeP.Q -0.021653 0.018686 -1.159 0.25157
Site.L:Genotype.L 0.033429 0.025724 1.300 0.19918
Site.L:Genotype.Q -0.007619 0.026324 -0.289 0.77333
Site.L:TimeP.L -0.039262 0.025724 -1.526 0.13267
Site.L:TimeP.Q -0.005353 0.026324 -0.203 0.83961
Genotype.L:TimeP.L 0.055269 0.031703 1.743 0.08686 .
Genotype.Q:TimeP.L 0.032487 0.031306 1.038 0.30393
Genotype.L:TimeP.Q -0.003858 0.031306 -0.123 0.90237
Genotype.Q:TimeP.Q -0.108144 0.033392 -3.239 0.00204 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.08798 on 55 degrees of freedom
Multiple R-squared: 0.7679, Adjusted R-squared: 0.7131
F-statistic: 14 on 13 and 55 DF, p-value: 4.691e-13
anova(Tol_ScoreF_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.18599 0.185989 24.0291 8.773e-06 ***
Genotype 2 0.58580 0.292900 37.8415 4.613e-11 ***
TimeP 2 0.48992 0.244959 31.6477 7.135e-10 ***
Site:Genotype 2 0.01497 0.007484 0.9669 0.38664
Site:TimeP 2 0.02075 0.010376 1.3405 0.27013
Genotype:TimeP 4 0.11112 0.027780 3.5891 0.01137 *
Residuals 55 0.42571 0.007740
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
----------------------------------------
Site | 0.10 | [0.01, 1.00]
Genotype | 0.32 | [0.15, 1.00]
TimeP | 0.27 | [0.10, 1.00]
Site:Genotype | 8.16e-03 | [0.00, 1.00]
Site:TimeP | 0.01 | [0.00, 1.00]
Genotype:TimeP | 0.06 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_lm.res<-data.frame(anova(Tol_ScoreF_lm))
Tol_ScoreF_lm.res$Predictor<-rownames(Tol_ScoreF_lm.res)
Tol_ScoreF_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_lm.res))
Tol_ScoreF_lm.res<-Tol_ScoreF_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Score_Full.prop)
shapiro.test(TolData_W2$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData_W2$Score_Full.prop
W = 0.96098, p-value = 0.4834
#Normal
##Model as a function of Site and Genotype
Tol_ScoreF_W2_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_W2_lm)); qqline(resid(Tol_ScoreF_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_W2_lm), resid(Tol_ScoreF_W2_lm))
#Model Results
summary(Tol_ScoreF_W2_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype,
data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.20185 -0.04952 -0.00045 0.04135 0.16725
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.706208 0.023330 30.270 3.15e-16 ***
Site.L 0.093421 0.032994 2.831 0.011517 *
Genotype.L -0.203806 0.040938 -4.978 0.000115 ***
Genotype.Q -0.068453 0.039874 -1.717 0.104193
Site.L:Genotype.L 0.003375 0.057895 0.058 0.954193
Site.L:Genotype.Q 0.010407 0.056391 0.185 0.855768
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1112 on 17 degrees of freedom
Multiple R-squared: 0.6882, Adjusted R-squared: 0.5964
F-statistic: 7.503 on 5 and 17 DF, p-value: 7e-04
anova(Tol_ScoreF_W2_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.12458 0.124577 10.0661 0.0055647 **
Genotype 2 0.33925 0.169625 13.7060 0.0002851 ***
Site:Genotype 2 0.00045 0.000226 0.0183 0.9819098
Residuals 17 0.21039 0.012376
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 0.18 | [0.00, 1.00]
Genotype | 0.50 | [0.17, 1.00]
Site:Genotype | 6.70e-04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_W2_lm.res<-data.frame(anova(Tol_ScoreF_W2_lm))
Tol_ScoreF_W2_lm.res$Predictor<-rownames(Tol_ScoreF_W2_lm.res)
Tol_ScoreF_W2_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_W2_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_W2_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_W2_lm.res))
Tol_ScoreF_W2_lm.res$TimeP<-rep("W2", nrow(Tol_ScoreF_W2_lm.res))
Tol_ScoreF_W2_lm.res<-Tol_ScoreF_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreF_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.755 0.0556 17 0.638 0.872
AC12 0.702 0.0556 17 0.585 0.819
AC8 0.463 0.0556 17 0.346 0.581
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.890 0.0556 17 0.772 1.007
AC12 0.822 0.0556 17 0.705 0.940
AC8 0.605 0.0642 17 0.469 0.740
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0529 0.0787 17 0.673 0.7820
AC10 - AC8 0.2916 0.0787 17 3.707 0.0047
AC12 - AC8 0.2387 0.0787 17 3.034 0.0195
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0676 0.0787 17 0.859 0.6723
AC10 - AC8 0.2848 0.0850 17 3.353 0.0100
AC12 - AC8 0.2172 0.0850 17 2.557 0.0508
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreF_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.755 0.0556 17 0.638 0.872
SS 0.890 0.0556 17 0.772 1.007
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.702 0.0556 17 0.585 0.819
SS 0.822 0.0556 17 0.705 0.940
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.463 0.0556 17 0.346 0.581
SS 0.605 0.0642 17 0.469 0.740
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.135 0.0787 17 -1.713 0.1049
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.120 0.0787 17 -1.527 0.1452
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.141 0.0850 17 -1.665 0.1142
##Save p-values
#Genotypes within Sites
Tol_ScoreF_W2_lm.geno<-data.frame(emmeans(Tol_ScoreF_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_W2_lm.geno<-Tol_ScoreF_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_W2_lm.geno$group1<-paste(Tol_ScoreF_W2_lm.geno$Site, Tol_ScoreF_W2_lm.geno$group1, sep="_")
Tol_ScoreF_W2_lm.geno$group2<-paste(Tol_ScoreF_W2_lm.geno$Site, Tol_ScoreF_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreF_W2_lm.site<-data.frame(emmeans(Tol_ScoreF_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_W2_lm.site<-Tol_ScoreF_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_W2_lm.site$group1<-paste(Tol_ScoreF_W2_lm.site$group1, Tol_ScoreF_W2_lm.site$Genotype, sep="_")
Tol_ScoreF_W2_lm.site$group2<-paste(Tol_ScoreF_W2_lm.site$group2, Tol_ScoreF_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreF_W2_lm.p<-rbind(Tol_ScoreF_W2_lm.geno[,c(1:2,4:8)], Tol_ScoreF_W2_lm.site[,c(1:2,4:8)])
Tol_ScoreF_W2_lm.p<-Tol_ScoreF_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreF_W2_lm.p$Sig<-ifelse(Tol_ScoreF_W2_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_W2_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreF_W2_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_W2_lm.p))
Tol_ScoreF_W2_lm.p$TimeP<-rep("W2", nrow(Tol_ScoreF_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreF_W2_SG<-summarySE(TolData_W2, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreF_W2_SG.plot<-ggplot(Tol_ScoreF_W2_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention Full Set")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_W2_lm.p, y.position=0.9, step.increase=0.35, label="Sig", hide.ns=TRUE); Tol_ScoreF_W2_SG.plot
##Check normality
hist(TolData_M1$Score_Full.prop)
shapiro.test(TolData_M1$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData_M1$Score_Full.prop
W = 0.89042, p-value = 0.01921
#Not normal
##Try square transformation
hist((TolData_M1$Score_Full.prop)^2)
shapiro.test((TolData_M1$Score_Full.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M1$Score_Full.prop)^2
W = 0.89427, p-value = 0.02287
#Not normal
##Try cubed transformation
hist((TolData_M1$Score_Full.prop)^3)
shapiro.test((TolData_M1$Score_Full.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M1$Score_Full.prop)^3
W = 0.88949, p-value = 0.01842
#Not normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreF_M1_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_M1_lm)); qqline(resid(Tol_ScoreF_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_M1_lm), resid(Tol_ScoreF_M1_lm))
#Model Results
summary(Tol_ScoreF_M1_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.205567 -0.027956 0.001925 0.027094 0.192833
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.83440 0.02043 40.846 < 2e-16 ***
Site.L 0.07153 0.02889 2.476 0.024843 *
Genotype.L -0.15896 0.03357 -4.736 0.000224 ***
Genotype.Q 0.08757 0.03711 2.360 0.031325 *
Site.L:Genotype.L 0.07642 0.04747 1.610 0.126961
Site.L:Genotype.Q 0.01094 0.05248 0.208 0.837562
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.09494 on 16 degrees of freedom
Multiple R-squared: 0.6982, Adjusted R-squared: 0.6039
F-statistic: 7.403 on 5 and 16 DF, p-value: 0.0009114
anova(Tol_ScoreF_M1_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.057569 0.057569 6.3867 0.0224118 *
Genotype 2 0.252333 0.126167 13.9969 0.0003061 ***
Site:Genotype 2 0.023755 0.011877 1.3177 0.2952999
Residuals 16 0.144222 0.009014
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.12 | [0.00, 1.00]
Genotype | 0.53 | [0.19, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_M1_lm.res<-data.frame(anova(Tol_ScoreF_M1_lm))
Tol_ScoreF_M1_lm.res$Predictor<-rownames(Tol_ScoreF_M1_lm.res)
Tol_ScoreF_M1_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_M1_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_M1_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M1_lm.res))
Tol_ScoreF_M1_lm.res$TimeP<-rep("M1", nrow(Tol_ScoreF_M1_lm.res))
Tol_ScoreF_M1_lm.res<-Tol_ScoreF_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreF_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.967 0.0475 16 0.866 1.068
AC12 0.719 0.0548 16 0.602 0.835
AC8 0.666 0.0475 16 0.565 0.766
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.998 0.0475 16 0.897 1.099
AC12 0.807 0.0548 16 0.691 0.923
AC8 0.850 0.0475 16 0.749 0.950
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.2484 0.0725 16 3.425 0.0092
AC10 - AC8 0.3012 0.0671 16 4.487 0.0010
AC12 - AC8 0.0528 0.0725 16 0.729 0.7504
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1909 0.0725 16 2.633 0.0452
AC10 - AC8 0.1484 0.0671 16 2.210 0.0998
AC12 - AC8 -0.0425 0.0725 16 -0.587 0.8292
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreF_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.967 0.0475 16 0.866 1.068
SS 0.998 0.0475 16 0.897 1.099
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.719 0.0548 16 0.602 0.835
SS 0.807 0.0548 16 0.691 0.923
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.666 0.0475 16 0.565 0.766
SS 0.850 0.0475 16 0.749 0.950
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0311 0.0671 16 -0.463 0.6499
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0885 0.0775 16 -1.142 0.2702
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.1839 0.0671 16 -2.739 0.0146
##Save p-values
#Genotypes within Sites
Tol_ScoreF_M1_lm.geno<-data.frame(emmeans(Tol_ScoreF_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_M1_lm.geno<-Tol_ScoreF_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M1_lm.geno$group1<-paste(Tol_ScoreF_M1_lm.geno$Site, Tol_ScoreF_M1_lm.geno$group1, sep="_")
Tol_ScoreF_M1_lm.geno$group2<-paste(Tol_ScoreF_M1_lm.geno$Site, Tol_ScoreF_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreF_M1_lm.site<-data.frame(emmeans(Tol_ScoreF_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_M1_lm.site<-Tol_ScoreF_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M1_lm.site$group1<-paste(Tol_ScoreF_M1_lm.site$group1, Tol_ScoreF_M1_lm.site$Genotype, sep="_")
Tol_ScoreF_M1_lm.site$group2<-paste(Tol_ScoreF_M1_lm.site$group2, Tol_ScoreF_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreF_M1_lm.p<-rbind(Tol_ScoreF_M1_lm.geno[,c(1:2,4:8)], Tol_ScoreF_M1_lm.site[,c(1:2,4:8)])
Tol_ScoreF_M1_lm.p<-Tol_ScoreF_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreF_M1_lm.p$Sig<-ifelse(Tol_ScoreF_M1_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_M1_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreF_M1_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M1_lm.p))
Tol_ScoreF_M1_lm.p$TimeP<-rep("M1", nrow(Tol_ScoreF_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreF_M1_SG<-summarySE(TolData_M1, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreF_M1_SG.plot<-ggplot(Tol_ScoreF_M1_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention Full Set")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_M1_lm.p, y.position=1.1, step.increase=0.2, label="Sig", hide.ns=TRUE); Tol_ScoreF_M1_SG.plot
##Check normality
hist(TolData_M4$Score_Full.prop)
shapiro.test(TolData_M4$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData_M4$Score_Full.prop
W = 0.83109, p-value = 0.0009936
#Not Normal
##Try square transformation
hist((TolData_M4$Score_Full.prop)^2)
shapiro.test((TolData_M4$Score_Full.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M4$Score_Full.prop)^2
W = 0.8317, p-value = 0.001019
#Not Normal
##Try cubed transformation
hist((TolData_M4$Score_Full.prop)^3)
shapiro.test((TolData_M4$Score_Full.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M4$Score_Full.prop)^3
W = 0.83161, p-value = 0.001015
#Not Normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreF_M4_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_M4_lm)); qqline(resid(Tol_ScoreF_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_M4_lm), resid(Tol_ScoreF_M4_lm))
#Model Results
summary(Tol_ScoreF_M4_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.137800 -0.019025 0.003675 0.020575 0.099500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.90906 0.01100 82.614 < 2e-16 ***
Site.L 0.03859 0.01556 2.480 0.0233 *
Genotype.L -0.12460 0.01906 -6.538 3.82e-06 ***
Genotype.Q -0.02191 0.01906 -1.149 0.2654
Site.L:Genotype.L 0.01901 0.02695 0.705 0.4896
Site.L:Genotype.Q -0.04168 0.02695 -1.546 0.1394
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05391 on 18 degrees of freedom
Multiple R-squared: 0.7468, Adjusted R-squared: 0.6765
F-statistic: 10.62 on 5 and 18 DF, p-value: 7.153e-05
anova(Tol_ScoreF_M4_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.017871 0.017871 6.1496 0.02326 *
Genotype 2 0.128043 0.064021 22.0308 1.452e-05 ***
Site:Genotype 2 0.008394 0.004197 1.4442 0.26199
Residuals 18 0.052308 0.002906
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.09 | [0.00, 1.00]
Genotype | 0.62 | [0.34, 1.00]
Site:Genotype | 0.04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_M4_lm.res<-data.frame(anova(Tol_ScoreF_M4_lm))
Tol_ScoreF_M4_lm.res$Predictor<-rownames(Tol_ScoreF_M4_lm.res)
Tol_ScoreF_M4_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_M4_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_M4_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M4_lm.res))
Tol_ScoreF_M4_lm.res$TimeP<-rep("M4", nrow(Tol_ScoreF_M4_lm.res))
Tol_ScoreF_M4_lm.res<-Tol_ScoreF_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreF_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.982 0.027 18 0.926 1.039
AC12 0.876 0.027 18 0.819 0.932
AC8 0.787 0.027 18 0.731 0.844
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.994 0.027 18 0.937 1.051
AC12 0.978 0.027 18 0.922 1.035
AC8 0.837 0.027 18 0.780 0.893
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1069 0.0381 18 2.804 0.0301
AC10 - AC8 0.1952 0.0381 18 5.122 0.0002
AC12 - AC8 0.0883 0.0381 18 2.318 0.0788
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0157 0.0381 18 0.411 0.9115
AC10 - AC8 0.1572 0.0381 18 4.124 0.0018
AC12 - AC8 0.1415 0.0381 18 3.713 0.0043
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreF_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.982 0.027 18 0.926 1.039
SS 0.994 0.027 18 0.937 1.051
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.876 0.027 18 0.819 0.932
SS 0.978 0.027 18 0.922 1.035
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.787 0.027 18 0.731 0.844
SS 0.837 0.027 18 0.780 0.893
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0115 0.0381 18 -0.302 0.7663
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.1027 0.0381 18 -2.694 0.0148
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.0495 0.0381 18 -1.299 0.2103
##Save p-values
#Genotypes within Sites
Tol_ScoreF_M4_lm.geno<-data.frame(emmeans(Tol_ScoreF_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_M4_lm.geno<-Tol_ScoreF_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M4_lm.geno$group1<-paste(Tol_ScoreF_M4_lm.geno$Site, Tol_ScoreF_M4_lm.geno$group1, sep="_")
Tol_ScoreF_M4_lm.geno$group2<-paste(Tol_ScoreF_M4_lm.geno$Site, Tol_ScoreF_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreF_M4_lm.site<-data.frame(emmeans(Tol_ScoreF_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_M4_lm.site<-Tol_ScoreF_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M4_lm.site$group1<-paste(Tol_ScoreF_M4_lm.site$group1, Tol_ScoreF_M4_lm.site$Genotype, sep="_")
Tol_ScoreF_M4_lm.site$group2<-paste(Tol_ScoreF_M4_lm.site$group2, Tol_ScoreF_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreF_M4_lm.p<-rbind(Tol_ScoreF_M4_lm.geno[,c(1:2,4:8)], Tol_ScoreF_M4_lm.site[,c(1:2,4:8)])
Tol_ScoreF_M4_lm.p<-Tol_ScoreF_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreF_M4_lm.p$Sig<-ifelse(Tol_ScoreF_M4_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_M4_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreF_M4_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M4_lm.p))
Tol_ScoreF_M4_lm.p$TimeP<-rep("M4", nrow(Tol_ScoreF_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreF_M4_SG<-summarySE(TolData_M4, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreF_M4_SG.plot<-ggplot(Tol_ScoreF_M4_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention Full Set")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_M4_lm.p, y.position=1.1, step.increase=0.22, label="Sig", hide.ns=TRUE); Tol_ScoreF_M4_SG.plot
##Check normality
hist(TolData$Score_TP.prop)
shapiro.test(TolData$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData$Score_TP.prop
W = 0.9161, p-value = 0.0001986
#Not Normal
##Try square transformation
hist((TolData$Score_TP.prop)^2)
shapiro.test((TolData$Score_TP.prop)^2)
Shapiro-Wilk normality test
data: (TolData$Score_TP.prop)^2
W = 0.92506, p-value = 0.0004895
#Not normal
##Try cubed transformation
hist((TolData$Score_TP.prop)^3)
shapiro.test((TolData$Score_TP.prop)^3)
Shapiro-Wilk normality test
data: (TolData$Score_TP.prop)^3
W = 0.92936, p-value = 0.0007674
#Not normal
##Model as a function of Site and Genotype and Timepoint
##Model with no transformation and check residuals
Tol_ScoreTP_lm<-lm(Score_TP.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_lm)); qqline(resid(Tol_ScoreTP_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_lm), resid(Tol_ScoreTP_lm))
Residuals are not great, need to check for other modeling options for Color by Timepoint.
#Model Results
summary(Tol_ScoreTP_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.133294 -0.019060 0.000725 0.026225 0.122906
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.883494 0.006319 139.809 < 2e-16 ***
Site.L 0.038546 0.008912 4.325 6.47e-05 ***
Genotype.L -0.102221 0.010797 -9.468 3.81e-13 ***
Genotype.Q 0.004971 0.011092 0.448 0.655784
TimeP.L 0.068495 0.010797 6.344 4.45e-08 ***
TimeP.Q 0.021265 0.011092 1.917 0.060411 .
Site.L:Genotype.L 0.014004 0.015269 0.917 0.363050
Site.L:Genotype.Q -0.008574 0.015625 -0.549 0.585424
Site.L:TimeP.L -0.036792 0.015269 -2.410 0.019344 *
Site.L:TimeP.Q -0.004886 0.015625 -0.313 0.755697
Genotype.L:TimeP.L 0.029281 0.018818 1.556 0.125448
Genotype.Q:TimeP.L 0.024642 0.018582 1.326 0.190292
Genotype.L:TimeP.Q 0.015340 0.018582 0.826 0.412653
Genotype.Q:TimeP.Q -0.072757 0.019821 -3.671 0.000548 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05222 on 55 degrees of freedom
Multiple R-squared: 0.763, Adjusted R-squared: 0.7069
F-statistic: 13.62 on 13 and 55 DF, p-value: 8.098e-13
anova(Tol_ScoreTP_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.057199 0.057199 20.9741 2.709e-05 ***
Genotype 2 0.238412 0.119206 43.7110 4.331e-12 ***
TimeP 2 0.117103 0.058552 21.4700 1.284e-07 ***
Site:Genotype 2 0.003221 0.001611 0.5906 0.557481
Site:TimeP 2 0.017283 0.008642 3.1687 0.049832 *
Genotype:TimeP 4 0.049580 0.012395 4.5450 0.003041 **
Residuals 55 0.149993 0.002727
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
----------------------------------------
Site | 0.09 | [0.01, 1.00]
Genotype | 0.38 | [0.20, 1.00]
TimeP | 0.19 | [0.04, 1.00]
Site:Genotype | 5.09e-03 | [0.00, 1.00]
Site:TimeP | 0.03 | [0.00, 1.00]
Genotype:TimeP | 0.08 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_lm.res<-data.frame(anova(Tol_ScoreTP_lm))
Tol_ScoreTP_lm.res$Predictor<-rownames(Tol_ScoreTP_lm.res)
Tol_ScoreTP_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_lm.res))
Tol_ScoreTP_lm.res<-Tol_ScoreTP_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with variables of interest, Site and Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Score_TP.prop)
shapiro.test(TolData_W2$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData_W2$Score_TP.prop
W = 0.95832, p-value = 0.4302
#Normal
##Model as a function of Site and Genotype
Tol_ScoreTP_W2_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_W2_lm)); qqline(resid(Tol_ScoreTP_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_W2_lm), resid(Tol_ScoreTP_W2_lm))
#Model Results
summary(Tol_ScoreTP_W2_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype,
data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.117075 -0.028442 -0.002075 0.023658 0.099425
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.843522 0.013080 64.488 < 2e-16 ***
Site.L 0.062257 0.018498 3.366 0.00367 **
Genotype.L -0.117129 0.022952 -5.103 8.83e-05 ***
Genotype.Q -0.042426 0.022356 -1.898 0.07485 .
Site.L:Genotype.L -0.002071 0.032459 -0.064 0.94987
Site.L:Genotype.Q 0.004044 0.031616 0.128 0.89972
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.06237 on 17 degrees of freedom
Multiple R-squared: 0.7169, Adjusted R-squared: 0.6336
F-statistic: 8.61 on 5 and 17 DF, p-value: 0.0003242
anova(Tol_ScoreTP_W2_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.054151 0.054151 13.9200 0.0016623 **
Genotype 2 0.113242 0.056621 14.5549 0.0002073 ***
Site:Genotype 2 0.000083 0.000041 0.0106 0.9894534
Residuals 17 0.066133 0.003890
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 0.23 | [0.01, 1.00]
Genotype | 0.48 | [0.15, 1.00]
Site:Genotype | 3.53e-04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_W2_lm.res<-data.frame(anova(Tol_ScoreTP_W2_lm))
Tol_ScoreTP_W2_lm.res$Predictor<-rownames(Tol_ScoreTP_W2_lm.res)
Tol_ScoreTP_W2_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_W2_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_W2_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_W2_lm.res))
Tol_ScoreTP_W2_lm.res$TimeP<-rep("W2", nrow(Tol_ScoreTP_W2_lm.res))
Tol_ScoreTP_W2_lm.res<-Tol_ScoreTP_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreTP_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.863 0.0312 17 0.797 0.929
AC12 0.836 0.0312 17 0.771 0.902
AC8 0.699 0.0312 17 0.633 0.765
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.955 0.0312 17 0.889 1.021
AC12 0.920 0.0312 17 0.854 0.986
AC8 0.788 0.0360 17 0.712 0.864
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0263 0.0441 17 0.597 0.8236
AC10 - AC8 0.1636 0.0441 17 3.709 0.0047
AC12 - AC8 0.1373 0.0441 17 3.112 0.0166
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0354 0.0441 17 0.803 0.7065
AC10 - AC8 0.1677 0.0476 17 3.521 0.0070
AC12 - AC8 0.1323 0.0476 17 2.778 0.0328
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreTP_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.863 0.0312 17 0.797 0.929
SS 0.955 0.0312 17 0.889 1.021
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.836 0.0312 17 0.771 0.902
SS 0.920 0.0312 17 0.854 0.986
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.699 0.0312 17 0.633 0.765
SS 0.788 0.0360 17 0.712 0.864
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0925 0.0441 17 -2.096 0.0513
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0834 0.0441 17 -1.890 0.0759
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.0883 0.0476 17 -1.854 0.0812
##Save p-values
#Genotypes within Sites
Tol_ScoreTP_W2_lm.geno<-data.frame(emmeans(Tol_ScoreTP_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_W2_lm.geno<-Tol_ScoreTP_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_W2_lm.geno$group1<-paste(Tol_ScoreTP_W2_lm.geno$Site, Tol_ScoreTP_W2_lm.geno$group1, sep="_")
Tol_ScoreTP_W2_lm.geno$group2<-paste(Tol_ScoreTP_W2_lm.geno$Site, Tol_ScoreTP_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreTP_W2_lm.site<-data.frame(emmeans(Tol_ScoreTP_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_W2_lm.site<-Tol_ScoreTP_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_W2_lm.site$group1<-paste(Tol_ScoreTP_W2_lm.site$group1, Tol_ScoreTP_W2_lm.site$Genotype, sep="_")
Tol_ScoreTP_W2_lm.site$group2<-paste(Tol_ScoreTP_W2_lm.site$group2, Tol_ScoreTP_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreTP_W2_lm.p<-rbind(Tol_ScoreTP_W2_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_W2_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_W2_lm.p<-Tol_ScoreTP_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreTP_W2_lm.p$Sig<-ifelse(Tol_ScoreTP_W2_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_W2_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreTP_W2_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_W2_lm.p))
Tol_ScoreTP_W2_lm.p$TimeP<-rep("W2", nrow(Tol_ScoreTP_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreTP_W2_SG<-summarySE(TolData_W2, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreTP_W2_SG.plot<-ggplot(Tol_ScoreTP_W2_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention by Timepoint")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_W2_lm.p, y.position=1, step.increase=0.45, label="Sig", hide.ns=TRUE); Tol_ScoreTP_W2_SG.plot
##Check normality
hist(TolData_M1$Score_TP.prop)
shapiro.test(TolData_M1$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData_M1$Score_TP.prop
W = 0.9028, p-value = 0.03383
#Not normal
##Try square transformation
hist((TolData_M1$Score_TP.prop)^2)
shapiro.test((TolData_M1$Score_TP.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M1$Score_TP.prop)^2
W = 0.90185, p-value = 0.03238
#Not normal
##Try cubed transformation
hist((TolData_M1$Score_TP.prop)^3)
shapiro.test((TolData_M1$Score_TP.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M1$Score_TP.prop)^3
W = 0.89818, p-value = 0.02735
#Not normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreTP_M1_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_M1_lm)); qqline(resid(Tol_ScoreTP_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_M1_lm), resid(Tol_ScoreTP_M1_lm))
#Model Results
summary(Tol_ScoreTP_M1_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.125200 -0.020725 -0.001675 0.015662 0.131000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.866131 0.012681 68.303 < 2e-16 ***
Site.L 0.041582 0.017933 2.319 0.034 *
Genotype.L -0.114746 0.020837 -5.507 4.78e-05 ***
Genotype.Q 0.064377 0.023036 2.795 0.013 *
Site.L:Genotype.L 0.050375 0.029467 1.710 0.107
Site.L:Genotype.Q 0.004277 0.032577 0.131 0.897
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05893 on 16 degrees of freedom
Multiple R-squared: 0.7446, Adjusted R-squared: 0.6647
F-statistic: 9.327 on 5 and 16 DF, p-value: 0.0002606
anova(Tol_ScoreTP_M1_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.019311 0.019311 5.5598 0.03144 *
Genotype 2 0.132460 0.066230 19.0683 5.822e-05 ***
Site:Genotype 2 0.010210 0.005105 1.4698 0.25940
Residuals 16 0.055573 0.003473
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.09 | [0.00, 1.00]
Genotype | 0.61 | [0.30, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_M1_lm.res<-data.frame(anova(Tol_ScoreTP_M1_lm))
Tol_ScoreTP_M1_lm.res$Predictor<-rownames(Tol_ScoreTP_M1_lm.res)
Tol_ScoreTP_M1_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_M1_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_M1_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M1_lm.res))
Tol_ScoreTP_M1_lm.res$TimeP<-rep("M1", nrow(Tol_ScoreTP_M1_lm.res))
Tol_ScoreTP_M1_lm.res<-Tol_ScoreTP_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreTP_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.968 0.0295 16 0.906 1.031
AC12 0.787 0.0340 16 0.715 0.859
AC8 0.755 0.0295 16 0.693 0.818
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.979 0.0295 16 0.917 1.041
AC12 0.841 0.0340 16 0.768 0.913
AC8 0.867 0.0295 16 0.805 0.930
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1815 0.0450 16 4.032 0.0026
AC10 - AC8 0.2127 0.0417 16 5.103 0.0003
AC12 - AC8 0.0312 0.0450 16 0.693 0.7710
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1385 0.0450 16 3.077 0.0187
AC10 - AC8 0.1119 0.0417 16 2.685 0.0408
AC12 - AC8 -0.0266 0.0450 16 -0.591 0.8269
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreTP_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.968 0.0295 16 0.906 1.031
SS 0.979 0.0295 16 0.917 1.041
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.787 0.0340 16 0.715 0.859
SS 0.841 0.0340 16 0.768 0.913
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.755 0.0295 16 0.693 0.818
SS 0.867 0.0295 16 0.805 0.930
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0109 0.0417 16 -0.262 0.7970
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0539 0.0481 16 -1.119 0.2795
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.1116 0.0417 16 -2.679 0.0165
##Save p-values
#Genotypes within Sites
Tol_ScoreTP_M1_lm.geno<-data.frame(emmeans(Tol_ScoreTP_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_M1_lm.geno<-Tol_ScoreTP_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M1_lm.geno$group1<-paste(Tol_ScoreTP_M1_lm.geno$Site, Tol_ScoreTP_M1_lm.geno$group1, sep="_")
Tol_ScoreTP_M1_lm.geno$group2<-paste(Tol_ScoreTP_M1_lm.geno$Site, Tol_ScoreTP_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreTP_M1_lm.site<-data.frame(emmeans(Tol_ScoreTP_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_M1_lm.site<-Tol_ScoreTP_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M1_lm.site$group1<-paste(Tol_ScoreTP_M1_lm.site$group1, Tol_ScoreTP_M1_lm.site$Genotype, sep="_")
Tol_ScoreTP_M1_lm.site$group2<-paste(Tol_ScoreTP_M1_lm.site$group2, Tol_ScoreTP_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreTP_M1_lm.p<-rbind(Tol_ScoreTP_M1_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_M1_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_M1_lm.p<-Tol_ScoreTP_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreTP_M1_lm.p$Sig<-ifelse(Tol_ScoreTP_M1_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_M1_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreTP_M1_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M1_lm.p))
Tol_ScoreTP_M1_lm.p$TimeP<-rep("M1", nrow(Tol_ScoreTP_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreTP_M1_SG<-summarySE(TolData_M1, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreTP_M1_SG.plot<-ggplot(Tol_ScoreTP_M1_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention by Timepoint")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_M1_lm.p, y.position=1.05, step.increase=0.3, label="Sig", hide.ns=TRUE); Tol_ScoreTP_M1_SG.plot
##Check normality
hist(TolData_M4$Score_TP.prop)
shapiro.test(TolData_M4$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData_M4$Score_TP.prop
W = 0.85397, p-value = 0.002593
#Not Normal
##Try square transformation
hist((TolData_M4$Score_TP.prop)^2)
shapiro.test((TolData_M4$Score_TP.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M4$Score_TP.prop)^2
W = 0.85439, p-value = 0.00264
#Not Normal
##Try cubed transformation
hist((TolData_M4$Score_TP.prop)^3)
shapiro.test((TolData_M4$Score_TP.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M4$Score_TP.prop)^3
W = 0.85457, p-value = 0.002661
#Not Normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreTP_M4_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_M4_lm)); qqline(resid(Tol_ScoreTP_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_M4_lm), resid(Tol_ScoreTP_M4_lm))
#Model Results
summary(Tol_ScoreTP_M4_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.070400 -0.009438 0.004100 0.011281 0.057100
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.940608 0.006243 150.654 < 2e-16 ***
Site.L 0.010536 0.008830 1.193 0.2483
Genotype.L -0.075254 0.010814 -6.959 1.68e-06 ***
Genotype.Q -0.007308 0.010814 -0.676 0.5078
Site.L:Genotype.L -0.006950 0.015293 -0.454 0.6549
Site.L:Genotype.Q -0.032086 0.015293 -2.098 0.0503 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03059 on 18 degrees of freedom
Multiple R-squared: 0.7531, Adjusted R-squared: 0.6846
F-statistic: 10.98 on 5 and 18 DF, p-value: 5.763e-05
anova(Tol_ScoreTP_M4_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.001332 0.0013321 1.4238 0.2483
Genotype 2 0.045732 0.0228662 24.4414 7.407e-06 ***
Site:Genotype 2 0.004311 0.0021557 2.3042 0.1285
Residuals 18 0.016840 0.0009355
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.02 | [0.00, 1.00]
Genotype | 0.67 | [0.41, 1.00]
Site:Genotype | 0.06 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_M4_lm.res<-data.frame(anova(Tol_ScoreTP_M4_lm))
Tol_ScoreTP_M4_lm.res$Predictor<-rownames(Tol_ScoreTP_M4_lm.res)
Tol_ScoreTP_M4_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_M4_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_M4_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M4_lm.res))
Tol_ScoreTP_M4_lm.res$TimeP<-rep("M4", nrow(Tol_ScoreTP_M4_lm.res))
Tol_ScoreTP_M4_lm.res<-Tol_ScoreTP_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreTP_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.989 0.0153 18 0.957 1.021
AC12 0.921 0.0153 18 0.888 0.953
AC8 0.890 0.0153 18 0.858 0.922
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.993 0.0153 18 0.960 1.025
AC12 0.973 0.0153 18 0.940 1.005
AC8 0.879 0.0153 18 0.847 0.911
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0686 0.0216 18 3.171 0.0139
AC10 - AC8 0.0995 0.0216 18 4.599 0.0006
AC12 - AC8 0.0309 0.0216 18 1.429 0.3478
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0199 0.0216 18 0.922 0.6335
AC10 - AC8 0.1134 0.0216 18 5.242 0.0002
AC12 - AC8 0.0934 0.0216 18 4.320 0.0011
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreTP_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.989 0.0153 18 0.957 1.021
SS 0.993 0.0153 18 0.960 1.025
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.921 0.0153 18 0.888 0.953
SS 0.973 0.0153 18 0.940 1.005
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.890 0.0153 18 0.858 0.922
SS 0.879 0.0153 18 0.847 0.911
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.00332 0.0216 18 -0.154 0.8795
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.05195 0.0216 18 -2.402 0.0273
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.01057 0.0216 18 0.489 0.6308
##Save p-values
#Genotypes within Sites
Tol_ScoreTP_M4_lm.geno<-data.frame(emmeans(Tol_ScoreTP_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_M4_lm.geno<-Tol_ScoreTP_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M4_lm.geno$group1<-paste(Tol_ScoreTP_M4_lm.geno$Site, Tol_ScoreTP_M4_lm.geno$group1, sep="_")
Tol_ScoreTP_M4_lm.geno$group2<-paste(Tol_ScoreTP_M4_lm.geno$Site, Tol_ScoreTP_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreTP_M4_lm.site<-data.frame(emmeans(Tol_ScoreTP_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_M4_lm.site<-Tol_ScoreTP_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M4_lm.site$group1<-paste(Tol_ScoreTP_M4_lm.site$group1, Tol_ScoreTP_M4_lm.site$Genotype, sep="_")
Tol_ScoreTP_M4_lm.site$group2<-paste(Tol_ScoreTP_M4_lm.site$group2, Tol_ScoreTP_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreTP_M4_lm.p<-rbind(Tol_ScoreTP_M4_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_M4_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_M4_lm.p<-Tol_ScoreTP_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreTP_M4_lm.p$Sig<-ifelse(Tol_ScoreTP_M4_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_M4_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreTP_M4_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M4_lm.p))
Tol_ScoreTP_M4_lm.p$TimeP<-rep("M4", nrow(Tol_ScoreTP_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreTP_M4_SG<-summarySE(TolData_M4, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreTP_M4_SG.plot<-ggplot(Tol_ScoreTP_M4_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention by Timepoint")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_M4_lm.p, y.position=1.1, step.increase=0.4, label="Sig", hide.ns=TRUE); Tol_ScoreTP_M4_SG.plot
Creating a heatmap to compare the direction and significance of pairwise comparisons across Tolerance metrics. Positive estimates indicate that Group 1 > Group 2. P values of less than 0.05 are considered significant.
##Combine Results
Tol_contrasts<-rbind(Tol_Chl_W2_lm.p, Tol_Chl_M1_lm.p, Tol_Chl_M4_lm.p,
Tol_Sym_W2_lm.p, Tol_Sym_M1_lm.p, Tol_Sym_M4_lm.p,
Tol_ScoreF_W2_lm.p, Tol_ScoreF_M1_lm.p, Tol_ScoreF_M4_lm.p,
Tol_ScoreTP_W2_lm.p, Tol_ScoreTP_M1_lm.p, Tol_ScoreTP_M4_lm.p)
##Create Contrast Column
Tol_contrasts$Contrast<-paste(Tol_contrasts$group1, Tol_contrasts$group2, sep=" v. ")
##Add Set column with Contrast and Timepoint
Tol_contrasts$Set<-paste(Tol_contrasts$TimeP, Tol_contrasts$Contrast, sep=" ")
##Re-order by Seasonal Timepoint
Tol_contrasts$TimeP<-factor(Tol_contrasts$TimeP, levels=c("W2", "M1", "M4"), ordered=TRUE)
Tol_contrasts<- Tol_contrasts %>% arrange(desc(as.numeric(TimeP)))
Tol_contrasts$Set<-factor(Tol_contrasts$Set, levels=c(unique(Tol_contrasts$Set)), ordered=TRUE)
Convert Estimate to the Response scale (instead of log +1 scale) for models where the response was log transformed
Tol_contrasts$Estimate<-Tol_contrasts$estimate
#Chl M4
Tol_contrasts$Estimate[which(Tol_contrasts$Response=="Chlorophyll" & Tol_contrasts$TimeP=="M4")]<-c(exp(Tol_contrasts$estimate[which(Tol_contrasts$Response=="Chlorophyll" & Tol_contrasts$TimeP=="M4")])-1)
#Sym M1
Tol_contrasts$Estimate[which(Tol_contrasts$Response=="Symbionts" & Tol_contrasts$TimeP=="M1")]<-c(exp(Tol_contrasts$estimate[which(Tol_contrasts$Response=="Symbionts" & Tol_contrasts$TimeP=="M1")])-1)
Tol_Pairs_W2.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="W2"),], aes(x=Response, y=Contrast, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts Summer 1")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_W2.plot
Tol_Pairs_M1.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="M1"),], aes(x=Response, y=Contrast, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts Summer 2")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_M1.plot
Tol_Pairs_M4.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="M4"),], aes(x=Response, y=Contrast, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts Month 4")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_M4.plot
Tol_Pairs.plot<-ggplot(data=Tol_contrasts, aes(x=Response, y=Set, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs.plot
Subsetting Contrasts to only include contrasts with significant differences in at least one of the thermal tolerance metrics (Chlorophyll or Symbiont or Color retention)
##Remove rows with non-significant p-values
Tol_contrasts_sig<-subset(Tol_contrasts, p < 0.05)
##Keep Sets where at least one metric shows a significant difference
Tol_Sets_sig<-data.frame(Set=c(unique(Tol_contrasts_sig$Set)))
Tol_contrasts_sig_Set<-merge(Tol_Sets_sig, Tol_contrasts, all.x=TRUE, all.y=FALSE)
Tol_Pairs_sig.plot<-ggplot(data=Tol_contrasts_sig_Set, aes(x=Response, y=Set, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Significant Pairwise Contrasts")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_sig.plot
Assessing the proportion of agreement between methods as the pairwise comparisons where two methods both find a significant difference out of the total pairwise comparisons. Calculating the proportion of agreement between Chlorophyll and Color Full Set, Chlorophyll and Color Timepoint, Symbionts and Color Full Set, Symbionts and Color Timepoint, and Chlorophyll and Symbionts. Comparing the proportion of agreement between methods of quantifying thermal tolerance.
names(Tol_contrasts)
[1] "group1" "group2" "estimate" "SE" "df" "t.ratio" "p" "Sig" "Response"
[10] "TimeP" "Contrast" "Set" "Estimate"
##Change long to short format
Tol_contrasts.agree<-Tol_contrasts %>% pivot_wider(names_from="Response", values_from=c("p", "Estimate"), id_cols=c("TimeP", "Contrast", "Set"))
##Convert p values to binary
#0: No significant difference (p>0.05) detected in that pairwise contrast
#1: Significant difference (p<0.05) detected in that pairwise contrast
Tol_contrasts.agree[,c(4:7)]<-ifelse(Tol_contrasts.agree[,c(4:7)]<0.05, 1, 0)
##Convert estimates to Greater or Less
#Greater: Positive estimates indicate A>B in the Contrast
#Less: Negative estimates indicate A<B in the Contrast
Tol_contrasts.agree[,c(8:11)]<-ifelse(Tol_contrasts.agree[,c(8:11)]<0, "Less", "Greater")
##Identify agreement
#1: Both metrics identified a significant difference in that pairwise contrast OR Neither metric identified a significant difference
#0: One metrics identified a significant difference in that pairwise contrast
#For each matching case (1), confirm matching Estimates
Tol_contrasts.agree$Chl.Sym<-ifelse(Tol_contrasts.agree$p_Chlorophyll==1 &
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Symbionts &
Tol_contrasts.agree$Estimate_Chlorophyll==Tol_contrasts.agree$Estimate_Symbionts, 1,
ifelse(Tol_contrasts.agree$p_Chlorophyll==0 &
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Symbionts, 1, 0))
Tol_contrasts.agree$Chl.ColF<-ifelse(Tol_contrasts.agree$p_Chlorophyll==1 & Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_FullSet &
Tol_contrasts.agree$Estimate_Chlorophyll==Tol_contrasts.agree$Estimate_Color_FullSet, 1,
ifelse(Tol_contrasts.agree$p_Chlorophyll==0 &
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_FullSet, 1, 0))
Tol_contrasts.agree$Chl.ColTP<-ifelse(Tol_contrasts.agree$p_Chlorophyll==1 &
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_TP & Tol_contrasts.agree$Estimate_Chlorophyll==Tol_contrasts.agree$Estimate_Color_TP, 1,
ifelse(Tol_contrasts.agree$p_Chlorophyll==0 &
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_TP, 1, 0))
Tol_contrasts.agree$Sym.ColF<-ifelse(Tol_contrasts.agree$p_Symbionts==1 &
Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_FullSet & Tol_contrasts.agree$Estimate_Symbionts==Tol_contrasts.agree$Estimate_Color_FullSet, 1,
ifelse(Tol_contrasts.agree$p_Symbionts==0 &
Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_FullSet, 1, 0))
Tol_contrasts.agree$Sym.ColTP<-ifelse(Tol_contrasts.agree$p_Symbionts==1 &
Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_TP & Tol_contrasts.agree$Estimate_Symbionts==Tol_contrasts.agree$Estimate_Color_TP, 1,
ifelse(Tol_contrasts.agree$p_Symbionts==0 &
Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_TP, 1, 0))
Tol_contrasts.agree$ColF.ColTP<-ifelse(Tol_contrasts.agree$p_Color_FullSet==1 &
Tol_contrasts.agree$p_Color_FullSet==Tol_contrasts.agree$p_Color_TP & Tol_contrasts.agree$Estimate_Color_FullSet==Tol_contrasts.agree$Estimate_Color_TP, 1,
ifelse(Tol_contrasts.agree$p_Color_FullSet==0 &
Tol_contrasts.agree$p_Color_FullSet==Tol_contrasts.agree$p_Color_TP, 1, 0))
##Convert to long to calculate instances of agreement
Tol_contrasts.agree.long<-Tol_contrasts.agree %>% pivot_longer(cols=c("Chl.Sym", "Chl.ColF", "Chl.ColTP", "Sym.ColF", "Sym.ColTP", "ColF.ColTP"), names_to="Compare", values_to="Agree")
##Sum Cases of Agreement
Tol_contrasts.agree.sum<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$Compare), sum)
names(Tol_contrasts.agree.sum)<-c("Compare", "Agree")
##Total Pairwise Comparisons
Tol_contrasts.agree.sum$Total<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$Compare), function(x) length(x))[,2]
prop.test(Tol_contrasts.agree.sum$Agree, Tol_contrasts.agree.sum$Total)
6-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.sum$Agree out of Tol_contrasts.agree.sum$Total
X-squared = 6.0509, df = 5, p-value = 0.3013
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3 prop 4 prop 5 prop 6
0.7777778 0.7777778 0.7037037 0.9259259 0.6666667 0.7407407
No significant difference in the proportion of agreement in detecting pairwise differences in thermal tolerance across comparisons of traditional bleaching or color score methods (Pearson’s chi-squared test p-value = 0.3013)
##Sum Cases of Agreement by Comparison and Timepoint
Tol_contrasts.agree.sum.tp<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$TimeP, Tol_contrasts.agree.long$Compare), sum)
names(Tol_contrasts.agree.sum.tp)<-c("TimeP", "Compare", "Agree")
##Total Pairwise Comparisons
Tol_contrasts.agree.sum.tp$Total<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$TimeP, Tol_contrasts.agree.long$Compare), function(x) length(x))[,3]
##Subset by Comparison
Tol_contrasts.agree.Chl.Sym<-subset(Tol_contrasts.agree.sum.tp, Compare=="Chl.Sym")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Chl.Sym$Agree, Tol_contrasts.agree.Chl.Sym$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.Chl.Sym$Agree out of Tol_contrasts.agree.Chl.Sym$Total
X-squared = 5.6842, df = 2, p-value = 0.0583
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
0.5555556 1.0000000 0.5555556
##Subset by Comparison
Tol_contrasts.agree.Chl.ColF<-subset(Tol_contrasts.agree.sum.tp, Compare=="Chl.ColF")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Chl.ColF$Agree, Tol_contrasts.agree.Chl.ColF$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.Chl.ColF$Agree out of Tol_contrasts.agree.Chl.ColF$Total
X-squared = 9, df = 2, p-value = 0.01111
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
1.0000000 0.8888889 0.4444444
##Subset by Comparison
Tol_contrasts.agree.Chl.ColTP<-subset(Tol_contrasts.agree.sum.tp, Compare=="Chl.ColTP")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Chl.ColTP$Agree, Tol_contrasts.agree.Chl.ColTP$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.Chl.ColTP$Agree out of Tol_contrasts.agree.Chl.ColTP$Total
X-squared = 9, df = 2, p-value = 0.01111
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
0.8888889 1.0000000 0.4444444
##Subset by Comparison
Tol_contrasts.agree.Sym.ColF<-subset(Tol_contrasts.agree.sum.tp, Compare=="Sym.ColF")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Sym.ColF$Agree, Tol_contrasts.agree.Sym.ColF$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.Sym.ColF$Agree out of Tol_contrasts.agree.Sym.ColF$Total
X-squared = 3, df = 2, p-value = 0.2231
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
0.5555556 0.8888889 0.5555556
##Subset by Comparison
Tol_contrasts.agree.Sym.ColTP<-subset(Tol_contrasts.agree.sum.tp, Compare=="Sym.ColTP")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Sym.ColTP$Agree, Tol_contrasts.agree.Sym.ColTP$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.Sym.ColTP$Agree out of Tol_contrasts.agree.Sym.ColTP$Total
X-squared = 5.0143, df = 2, p-value = 0.0815
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
0.6666667 1.0000000 0.5555556
##Subset by Comparison
Tol_contrasts.agree.ColF.ColTP<-subset(Tol_contrasts.agree.sum.tp, Compare=="ColF.ColTP")
##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.ColF.ColTP$Agree, Tol_contrasts.agree.ColF.ColTP$Total)
Warning: Chi-squared approximation may be incorrect
3-sample test for equality of proportions without continuity correction
data: Tol_contrasts.agree.ColF.ColTP$Agree out of Tol_contrasts.agree.ColF.ColTP$Total
X-squared = 1.08, df = 2, p-value = 0.5827
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3
0.8888889 0.8888889 1.0000000
Tol_Chl_W2_SG.plot.fig<-ggplot(Tol_Chl_W2_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Summer 1")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="top",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Chlorophyll Retained")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_W2_lm.p, y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE)
Tol_Chl_M1_SG.plot.fig<-ggplot(Tol_Chl_M1_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Summer 2")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="top",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_M1_lm.p, y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE)
Tol_Chl_M4_SG.plot.fig<-ggplot(Tol_Chl_M4_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Winter")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="top",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_M4_lm.p, y.position=0.55, step.increase=0.20, label="Sig", hide.ns=TRUE)
Tol_Sym_W2_SG.plot.fig<-ggplot(Tol_Sym_W2_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.position="none")+
labs(x="Site and Genotype", y="Proportion Symbionts Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_W2_lm.p, y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE)
Tol_Sym_M1_SG.plot.fig<-ggplot(Tol_Sym_M1_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.position="none")+
labs(x="Site and Genotype", y="")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_M1_lm.p, y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE)
Tol_Sym_M4_SG.plot.fig<-ggplot(Tol_Sym_M4_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.position="none")+
labs(x="Site and Genotype", y="")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_M4_lm.p, y.position=1.1, step.increase=0.15, label="Sig", hide.ns=TRUE)
Tol_ScoreF_W2_SG.plot.fig<-ggplot(Tol_ScoreF_W2_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.position="none")+
labs(x="", y="Proportion Color Retained (Full)")+
ylim(0.25, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_W2_lm.p, y.position=0.9, step.increase=0.35, label="Sig", hide.ns=TRUE)
Tol_ScoreF_M1_SG.plot.fig<-ggplot(Tol_ScoreF_M1_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.position="none")+
labs(x="", y="")+
ylim(0.25, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_M1_lm.p, y.position=1.1, step.increase=0.2, label="Sig", hide.ns=TRUE)
Tol_ScoreF_M4_SG.plot.fig<-ggplot(Tol_ScoreF_M4_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.position="none")+
labs(x="", y="")+
ylim(0.25, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_M4_lm.p, y.position=1.1, step.increase=0.22, label="Sig", hide.ns=TRUE)
Tol_ScoreTP_W2_SG.plot.fig<-ggplot(Tol_ScoreTP_W2_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.position="none")+
labs(x="", y="Proportion Color Retained (Split)")+
ylim(0.25, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_W2_lm.p, y.position=1, step.increase=0.45, label="Sig", hide.ns=TRUE)
Tol_ScoreTP_M1_SG.plot.fig<-ggplot(Tol_ScoreTP_M1_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.position="none")+
labs(x="", y="")+
ylim(0.25, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_M1_lm.p, y.position=1.05, step.increase=0.3, label="Sig", hide.ns=TRUE)
Tol_ScoreTP_M4_SG.plot.fig<-ggplot(Tol_ScoreTP_M4_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.position="none")+
labs(x="", y="")+
ylim(0.25, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_M4_lm.p, y.position=1.1, step.increase=0.4, label="Sig", hide.ns=TRUE)
##Create Panel
Tolerance_fig<-plot_grid(Tol_Chl_W2_SG.plot.fig, Tol_Chl_M1_SG.plot.fig, Tol_Chl_M4_SG.plot.fig,
Tol_ScoreF_W2_SG.plot.fig, Tol_ScoreF_M1_SG.plot.fig, Tol_ScoreF_M4_SG.plot.fig,
Tol_ScoreTP_W2_SG.plot.fig, Tol_ScoreTP_M1_SG.plot.fig, Tol_ScoreTP_M4_SG.plot.fig,
Tol_Sym_W2_SG.plot.fig, Tol_Sym_M1_SG.plot.fig, Tol_Sym_M4_SG.plot.fig,
nrow=4, ncol=3, byrow=TRUE,
rel_widths=1,
rel_heights=c(1, 0.85, 0.85, 0.85),
labels=c("A", "B", "C",
"D", "E", "F",
"G", "H", "I",
"J", "K", "L"),
label_size=panel.lab.sz,
label_fontface = "bold")
##Save Figure
ggsave(filename="Figures/FigureS2_Tolerance_Across_Metrics.png", plot=Tolerance_fig, dpi=300, width=14, height=20, units="in")
##Save Figure
ggsave(filename="Figures/Figure4_Tolerance_Heatmap.png", plot=Tol_Pairs.plot, dpi=300, width=8, height=12, units="in")
##Combine Results Tables
TableS2_Tol.lm.res<-rbind(Tol_Chl_W2_lm.res, Tol_Chl_M1_lm.res, Tol_Chl_M4_lm.res,
Tol_Sym_W2_lm.res, Tol_Sym_M1_lm.res, Tol_Sym_M4_lm.res,
Tol_ScoreF_W2_lm.res, Tol_ScoreF_M1_lm.res, Tol_ScoreF_M4_lm.res,
Tol_ScoreTP_W2_lm.res, Tol_ScoreTP_M1_lm.res, Tol_ScoreTP_M4_lm.res)
##Add a Season Variable
TableS2_Tol.lm.res$Season<-ifelse(TableS2_Tol.lm.res$TimeP=="W2", "Summer1", ifelse(TableS2_Tol.lm.res$TimeP=="M1", "Summer2", ifelse(TableS2_Tol.lm.res$TimeP=="M4", "Winter" , NA)))
##Organize
names(TableS2_Tol.lm.res)
TableS2_Tol.lm.res<-TableS2_Tol.lm.res[,c("TimeP", "Season", "Response", "Predictor", "DF", "Sum.Sq", "Mean.Sq", "F.value", "p.value", "EtaSq")]
#Round to 4 digits
TableS2_Tol.lm.res$Sum.Sq<-round(TableS2_Tol.lm.res$Sum.Sq, 4)
TableS2_Tol.lm.res$Mean.Sq<-round(TableS2_Tol.lm.res$Mean.Sq, 4)
TableS2_Tol.lm.res$F.value<-round(TableS2_Tol.lm.res$F.value, 4)
TableS2_Tol.lm.res$p.value<-round(TableS2_Tol.lm.res$p.value, 4)
TableS2_Tol.lm.res$EtaSq<-round(TableS2_Tol.lm.res$EtaSq, 4)
##Write Out Table
write.csv(TableS2_Tol.lm.res, "Tables/TableS2_Tolerance_LM_Results.csv", row.names=FALSE)
##Pairwise Results Table
TableS3_Tol.pairwise<-Tol_contrasts
##Add a Season Variable
TableS3_Tol.pairwise$Season<-ifelse(TableS3_Tol.pairwise$TimeP=="W2", "Summer1", ifelse(TableS3_Tol.pairwise$TimeP=="M1", "Summer2", ifelse(TableS3_Tol.pairwise$TimeP=="M4", "Winter" , NA)))
##Re-order by Seasonal Timepoint
TableS3_Tol.pairwise<- TableS3_Tol.pairwise %>% arrange(as.numeric(TimeP))
##Organize
names(TableS3_Tol.pairwise)
TableS3_Tol.pairwise<-TableS3_Tol.pairwise[,c("TimeP", "Season", "Response", "Contrast", "Estimate", "SE", "df", "t.ratio", "p")]
TableS3_Tol.pairwise<-TableS3_Tol.pairwise %>% dplyr::rename( DF = df) %>% dplyr::rename( p.value = p)
#Round to 4 digits
TableS3_Tol.pairwise$Estimate<-round(TableS3_Tol.pairwise$Estimate, 4)
TableS3_Tol.pairwise$SE<-round(TableS3_Tol.pairwise$SE, 4)
TableS3_Tol.pairwise$t.ratio<-round(TableS3_Tol.pairwise$t.ratio, 4)
TableS3_Tol.pairwise$p.value<-round(TableS3_Tol.pairwise$p.value, 4)
##Write Out Table
write.csv(TableS3_Tol.pairwise, "Tables/TableS3_Tolerance_Pairwise_Results.csv", row.names=FALSE)